#######################################
#########FACTORIAL EXPERIMENT##########
#######################################

#LOAD /INSTALL PACKAGES
pacman::p_load(
  ggplot2, aod, gridExtra, cregg, margins, prediction, foreign, binom, Barnard, haven, Rmisc, sjmisc, magrittr, tidyverse, jtools, MASS,
  ggsci, 
  patchwork
)

####ANALYSIS WITH SPEEDERS####

### SET WORKING DIRECTORY TO FILE LOCATION
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#LOAD DATA
paper_risp_f<-read.dta("immigration_speeders.dta")

#EXPERIMENTAL FACTORS
paper_risp_f$approved_asy  <- factor(paper_risp_f$approved_asy)
paper_risp_f$exp_1  <- factor(paper_risp_f$exp_1)
paper_risp_f$exp_2  <- factor(paper_risp_f$exp_2)
paper_risp_f$exp_3  <- factor(paper_risp_f$exp_3)
paper_risp_f$ideo  <- factor(paper_risp_f$ideo)
paper_risp_f$ideo2  <- factor(paper_risp_f$ideo2)


#GLM 3WAY TREATMENT INTERACTIONS (NOT WEIGHTED)
m3way <- glm(approved_asy ~ exp_1*exp_2*exp_3, family=binomial(), data = paper_risp_f)
  summary(m3way)

##PREDICTIONS 3WAY
prediction(m3way, at = list(exp_2 = c("Looking for job", "Fleeing from war"), exp_1 = c("Low skilled", "Skilled"), exp_3 = c("Syrian", "Ukrainian"))) %>% 
  summary -> pred3way

### GGPLOT 3WAY
pred3way %>% 
  #filter(`at(exp_3)` == "Ukrainian") %>%
  ggplot(data = .,
         aes(x = `at(exp_1)`,
             y = Prediction,
             shape = `at(exp_2)`,
             color = `at(exp_2)`)) +
  facet_grid(. ~ `at(exp_3)`) + 
  geom_point(size = 4, position = position_dodge(width = 0)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0) # 
  ) +
  #ggsci::scale_color_lancet(name = "Description of at(exp_1)") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("firebrick2", "cornflowerblue")) +
  scale_shape_manual(name="", values = c(16, 17)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18)) +
  theme(legend.title = element_text(size = 18),
    legend.text = element_text(size = 18)
  ) +
  theme(strip.text = element_text(size=18)) 
  

ggsave("figure3.png", dpi = 600)



#GLM 2  MAIN EFFECTS BY IDEOLOGY (NOT WEIGHTED)
m2way <- glm(approved_asy ~ exp_1*ideo2+exp_2*ideo2+exp_3*ideo2, family=binomial(), data = paper_risp_f)
summary(m2way)

#PREDICTIONS 2WAY
prediction(m2way, at = list(ideo2 = c("left", "right"), exp_1 = c("Low skilled", "Skilled"))) %>% 
  summary -> pred2wayskills

prediction(m2way, at = list(ideo2 = c("left", "right"), exp_2 = c("Looking for job", "Fleeing from war"))) %>% 
  summary -> pred2wayreason

prediction(m2way, at = list(ideo2 = c("left", "right"), exp_3 = c("Syrian", "Ukrainian"))) %>% 
  summary -> pred2wayreasonorigin

### GGPLOT 2WAY
ggplot(data = pred2wayskills,
       aes(x = `at(exp_1)`,
           y = Prediction,
           shape = `at(ideo2)`,
           color = `at(ideo2)`)) +
  geom_point(size = 4, position = position_dodge(width = 0.5)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0.5)
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("firebrick2", "cornflowerblue")) +
  scale_shape_manual(name="", values = c(16, 17)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18,face="bold")) +
  theme(legend.title = element_text(size = 20),
        legend.text = element_text(size = 20)
  ) +
  theme(strip.text = element_text(size=18)) -> plot1

ggplot(data = pred2wayreason,
       aes(x = `at(exp_2)`,
           y = Prediction,
           shape = `at(ideo2)`,
           color = `at(ideo2)`)) +
  geom_point(size = 4, position = position_dodge(width = 0.5)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0.5)
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("firebrick2", "cornflowerblue")) +
  scale_shape_manual(name="", values = c(16, 17)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18,face="bold")) +
  theme(legend.title = element_text(size = 20),
        legend.text = element_text(20)
  ) +
  theme(strip.text = element_text(size=18)) -> plot2

ggplot(data = pred2wayreasonorigin,
       aes(x = `at(exp_3)`,
           y = Prediction,
           shape = `at(ideo2)`,
           color = `at(ideo2)`)) +
  geom_point(size = 4, position = position_dodge(width = 0.5)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0.5)
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("firebrick2", "cornflowerblue")) +
  scale_shape_manual(name="", values = c(16, 17)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18,face="bold")) +
  theme(legend.title = element_text(size = 20),
        legend.text = element_text(size = 20)
  ) +
  theme(strip.text = element_text(size=18)) -> plot3


g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

mylegend<-g_legend(plot1)
plotcombined <- grid.arrange(arrangeGrob(plot1 + theme(legend.position="none"),
                                         plot2 + theme(legend.position="none"),
                                         plot3 + theme(legend.position="none"),
                                         nrow=1),
                             mylegend, nrow=2,heights=c(10, 1))

ggsave("figure4.png", dpi = 600, plotcombined)
